# Set working directory and log file
setwd("/Users/fredsolt/Documents/Projects/Religiosity/ Chapter")

# Include required packages
library(foreign)

# Enter data
ch2_res <- read.dta("ch2_res.dta", convert.factors=FALSE, convert.underscore=FALSE) 
gini_res <- read.dta("gini_net_res.dta", convert.factors=FALSE, convert.underscore=FALSE)
rgdpl_res <- read.dta("rgdpl_res.dta", convert.factors=FALSE, convert.underscore=FALSE) 

#####Granger-Test Graph

##point estimates, in a n.variables, n.variables x n.models
coef.matrix <- matrix(c(
						ch2_res$coef_all,
						gini_res$coef_all,
						rgdpl_res$coef_all
						),nr=4)

##standard error matrix, n.variables x n.models
se.matrix <- matrix(c(
						ch2_res$se_all,
						gini_res$se_all,
						rgdpl_res$se_all
                      ),nr=4)

##variable names
varnames<- c("Lagged Inequality","Lagged GDP/Capita","Lagged Religiosity", "Constant")

##exclude intercept
coef.matrix<-coef.matrix[-4,]
se.matrix<-se.matrix[-4,]

##we are making a list, define it first as empty
Y1 <- vector(length=0,mode="list")
#estimates
Y1$estimate <- coef.matrix
##95% confidence intervals
Y1$lo <- coef.matrix-qnorm(0.975)*se.matrix
Y1$hi <- coef.matrix+qnorm(0.975)*se.matrix
##90% confidence intervals
Y1$lo1 <- coef.matrix-qnorm(0.95)*se.matrix
Y1$hi1 <- coef.matrix+qnorm(0.95)*se.matrix

# ##we are making a list, define it first as empty
# Y1 <- vector(length=0,mode="list")
# #estimates
# Y1$estimate <- coef.matrix[1:2,]
# ##95% confidence intervals
# Y1$lo <- coef.matrix[1:2,]-qnorm(0.975)*se.matrix[1:2,]
# Y1$hi <- coef.matrix[1:2,]+qnorm(0.975)*se.matrix[1:2,]
# ##90% confidence intervals
# Y1$lo1 <- coef.matrix[1:2,]-qnorm(0.95)*se.matrix[1:2,]
# Y1$hi1 <- coef.matrix[1:2,]+qnorm(0.95)*se.matrix[1:2,]


##name the rows of Y1 estimate
rownames(Y1$estimate) <- varnames[-4]

##Plot estimates in a single columns
source("plotReg.R")
pdf(file="coef_graph.pdf",width=5.5,height=3)
## create the graph (do not print it yet)
tmp <- plot.reg(Y1,#the lists
                #the model labels
                label.x=c("Religiosity","Inequality","GDP/Capita"),
                ## reference lines
                refline=rep(0, times=3),
                ## space left in the bottom (for the x-axis labels)
                hlast=.2,
                ## print the graph?
                print=F,
                ## line width / character size multiplier
                lwd.fact=1,
                ## length of the cross-hairs
                length.arrow=unit(0,"mm"),
                ## widths: variable names, plot size, y-axis
                widths=c(.25,.65,.1),
                ## rotation of the variable name labels
                rot.label.y=0,
                ## justification of the variable name labels
                just.label.y="right",
                ## position (x-axis) of the variable name labels)
                pos.label.y=0.95,
                ## size of the symbol
                pch.size=0.2,expand.factor=.2,expand.factor2=0.1,
                v.grid=T,
                ## y-axes
                yaxis.at=list(
                seq(0,1,.5),
                seq(-1,1,1),
                seq(0,1,.5)
                )
                )
grid.draw(tmp) ## print the graph
## Additional Info
# fit <- textGrob(expression(R^2), x=.23, y=.23, just="right", gp=gpar(cex=.8))
# fit1 <- textGrob(".251", x=.355, y=.23, just="center", gp=gpar(cex=.8))
# fit2 <- textGrob(".832", x=.57, y=.23, just="center", gp=gpar(cex=.8))
# fit3 <- textGrob(".985", x=.79, y=.23, just="center", gp=gpar(cex=.8))
# grid.draw(fit)
# grid.draw(fit1)
# grid.draw(fit2)
# grid.draw(fit3)
# fit.b <- textGrob(expression(R^2), x=.23, y=.18, just="right", gp=gpar(cex=.77))
# fit.b1 <- textGrob(".847", x=.355, y=.18, just="center", gp=gpar(cex=.8))
# fit.b2 <- textGrob(".998", x=.57, y=.18, just="center", gp=gpar(cex=.8))
# fit.b3 <- textGrob(".999", x=.79, y=.18, just="center", gp=gpar(cex=.8))
# grid.draw(fit.b)
# grid.draw(fit.b1)
# grid.draw(fit.b2)
# grid.draw(fit.b3)
fit.o <- textGrob(expression(R^2), x=.23, y=.07, just="right", gp=gpar(cex=.8))
fit.o1 <- textGrob(".780", x=.355, y=.07, just="center", gp=gpar(cex=.8))
fit.o2 <- textGrob(".977", x=.57, y=.07, just="center", gp=gpar(cex=.8))
fit.o3 <- textGrob(".995", x=.79, y=.07, just="center", gp=gpar(cex=.8))
grid.draw(fit.o)
grid.draw(fit.o1)
grid.draw(fit.o2)
grid.draw(fit.o3)
n <- textGrob("N", x=.23, y=.18, just="right", gp=gpar(cex=.8))
n1 <- textGrob("788", x=.355, y=.18, just="center", gp=gpar(cex=.8))
n2 <- textGrob("787", x=.57, y=.18, just="center", gp=gpar(cex=.8))
n3 <- textGrob("787", x=.79, y=.18, just="center", gp=gpar(cex=.8))
grid.draw(n)
grid.draw(n1)
grid.draw(n2)
grid.draw(n3)
graphics.off()

q()

